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Abstract 

We describe the construction of a fully tractable mathematical model for intracellular pH. This work is based on coupling the 
kinetic equations depicting the molecular mechanisms for pumps, transporters and chemical reactions, which determine 
this parameter in eukaryotic cells. Thus, our system also calculates the membrane potential and the cytosolic ionic 
composition. Such a model required the development of a novel algebraic method that couples differential equations for 
slow relaxation processes to steady-state equations for fast chemical reactions. Compared to classical heuristic approaches 
based on fitted curves and ad hoc constants, this yields significant improvements. This model is mathematically self- 
consistent and allows for the first time to establish analytical solutions for steady-state pH and a reduced differential 
equation for pH regulation. Because of its modular structure, it can integrate any additional mechanism that will directly or 
indirectly affect pH. In addition, it provides mathematical clarifications for widely observed biological phenomena such as 
overshooting in regulatory loops. Finally, instead of including a limited set of experimental results to fit our model, we show 
examples of numerical calculations that are extremely consistent with the wide body of intracellular pH experimental 
measurements gathered by different groups in many different cellular systems. 
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Introduction 

Distribution of charges within biological molecules is crucial, 
not only for reactivity and catalysis, but also as it determines their 
solubility, their particular folding, and dictates the spatio-temporal 
sequence of their interactions. In this context, the pH of the 
solution bathing these biological molecules is a key parameter, 
since its value determines the protonation of the acid-base groups 
that are especially abundant in macromolecular assemblies. 
Furthermore, as many enzymes and cellular regulators exhibit a 
strong pH dependency, the modification of the protonation of key 
residues can deeply impact their function. For these reasons, 
genomes necessarily contain pH-dependency information, which is 
expressed in the proteome [1]. The complete information for 
intracellular pH determination is a convoluted interplay between 
the abundance and the distribution of protonable groups in 
biological molecules, their pKa values and the expression, stability, 
kinetic and affinity parameters of the pH regulating systems. 
Accordingly, providing a fully tractable model for intracellular pH 
regulation is a challenging problem, and several studies have been 
aimed at building essentially heuristic models [2-5] for intracel- 
lular pH regulation. 

The past decades have witnessed the detailed molecular 
characterization of the protagonists that regulate the concentra- 
tions of cellular acid-base equivalents, in term of both their kinetics 



and the affinities for their substrates [6,7]. Significant efforts have 
also been invested to describe intracellular buffering mechanisms 
and proton diffusion in cells adequately [8,9]. 

Based on this, we develop here a different, bottom-up approach 
at the interface between biology, physics, chemistry and mathe- 
matics. We construct a model that encompasses the individual 
molecular mechanisms for these regulators defined by their own 
kinetics and by their experimentally measured microscopic 
parameters. This requires the inclusion of the chemical reactions 
between the involved reactive species. This non-empirical process 
guarantees the construction of a physically coherent, fully 
integrated and tractable model (i) for cellular proton dynamics 
and (ii) for steady-state pH regulation. 

In the present study, we choose to keep the system simple and 
modular by assuming that the cell surface and volume are fixed to 
their average values and by using the ubiquitous Na + /H + 
exchanger NHE— 1 and Cl^/HCO^ exchanger AE2 as the 
main transmembrane acid-base transporters. We also include the 
electrical gradient generated by the Na/K-ATPase across the 
membrane and the permeabilities associated to Na + , K + and 
CI - background currents measured in non-excitable cells. 
Therefore, our model computes the distribution of the other 
cationic and anionic species and their variations as a function of 
proton concentration. 
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These pumps and transporters show a very high sequence 
conservation within different mammalian species and possess very 
similar constants for their substrates. Based on this, we built our 
model using widely accepted values from the literature even if they 
had been measured from different mammalian species. We will 
further see that this is validated by our results, which show that pH 
regulation is very resilient against variations of those thermody- 
namic constants. 

It is demonstrated that our model gives (i) a robust, experiment 
based prediction of the temporal evolution of the pH, (ii) a simple 
analytical value for its steady state, (hi) all the other ionic 
concentrations related to the proton regulation, (iv) and a reduced 
differential equation for describing the full pH balance. 

This enables the testing of biologically-relevant situations whilst 
discriminating between critical parameters and rate limiting steps 
versus those factors that can be widely changed with virtually no 
effect on cellular homeostasis. 

Methods 

Datasets used for the Model 

We report most of thermodynamical data, the common ionic 
environments, and the justification of the kinetic equations in the 
Datasets in File SI. In the following, we illustrate the specific 
behavior of the involved physical, chemical or biological 
components. 

Ionic Flows and Potential through the Membrane 

Let us depict the cellular model represented in Figure 1 
mathematically. We assume that the cell geometry is fixed by 
neglecting that water flows through the membrane. The charge 
balance is controlled by passive, electroneutral, electrogenic flows 
and capacitive currents that are described as follows. 

Passive flows. If X represents a chemical species in Figure 1, 
with an inner concentration [X] and an outer concentration [X] out , 
then it flows out of the cell through the membrane surface S due to 
a permeability Px(£m)- Here, C m = J-E„,/(RT) represents the 
reduced electric potential where T is the Faraday constant, R is 
the molar gas constant, T is the absolute temperature, and E m is 
the electric potential difference between the cytosol and the outer 
medium. The Goldman-Hodgkin-Katz flux equation [10] pro- 
vides the outward molar flux jx as 



-iMU^xQ ([x] out - [X]^ 



with 'F(m) = u{e" — 1) and zx is the algebraic charge of X. The 
associated passive outward electric flux is Jx = z x3~jx and the 
whole cell passive outward electric current is Ix = SJx- We can 
simply convert the flux into an intake molar rate for a given cell 
volume V as 



3, [XI 



V 



7x = ^x(C m ^x,[X] out ,[X]). 



(2) 



For the cellular system used in the electrophysiological measure- 
ments we recorded significant currents only for K + , Na + and 
Cl~ (CCL39 cells, see Figure SI in File SI). This allows the 
determination of the corresponding permeabilities. Any other 
species can be taken into account if other cells are considered, and 
if values are available or measurable. 

Electroneutral transports. The electroneutral AE2 ex- 
changer keeps Cl~ ion concentration above its Nernst potential 



[1 1], and thus is assumed to work in the forward direction. AE2 is 
then operating with a Hill mechanism [12], inducing a whole cell 
exchange rate 



3 f [Cri =-5 ( [HC0 3 



= Pae 



(3) 



with Pae =k ae [hco 3 -]7g4 e 

the cellular maximal HC0 3 /CI" 



[HCO3-] 1 ') and where V AE is 
exchange rate, and where -Kae 
is the bicarbonate affinity. Unless indicated, we will use a 
Michaelis-Menten behavior, namely y = 1 , and -Kae is about 
10 mM. 

We use the established mechanism [13] of the NHE— 1 
exchanger that results in the whole cell exchange rate 



5 <[ Na+ UE=-^ H+ U E ^NHE 



(4) 



with p NH E= ^ / nheO'nhe([H + ]/A:,.) and where K NH e is the 
cellular maximum H + /Na + exchange rate, K r is about 
1.810 8 M and 



cnhe(-v-) : 



x(\+x) + L 0 Cx(l + Cx) 
L 0 (l + Cx) 2 + (l+x) 2 



with C = K r /K h K,~3.6-\0- 6 M and L 0 ^10 3 . 

Any other electroneutral transporter could be similarly 
described and therefore inserted into the model. 

Electrogenic currents. We restrict ourselves to the sodium- 
potassium pump that exchanges three inner Na + with two outer 
K + according to 



^,[K+[ 

2 L |NaK 



3 5 f[Na ] |NaK =PNaK&«). 



(5) 



where p NaK (Cm) = ^NaKffNaK(C m ,[Na + [). K NaK is the cellular 
maximum Na + /K + exchange rate, and we combined the 
experimental data found in published studies [14,15] to estimate 
(see Datasets in File SI) 



ffNaK(C„„[Na+]) = 



[1+ tanh(0.39C,„+1.28)[ [Na~ 



^NaK + [Na + ] 



with ^T NaK «10— 15 mM. 

Electric potential evolution. The cytosol and the outer 
medium must remain globally electroneutral. Conversely, charge 
accumulation polarizes the membrane due to its surface capac- 
itance C m ^ 1 ^F.cm~ 2 (see Materials and Methods in File SI). 
The total capacitance of cell membrane is C = SC,„ . 

We take into account both the passive and electrogenic actors, 
respectively defined by equations (1) and (5), which are involved in 
the electric potential regulation. This results in 



Cd,E m + 7 x(£,„) + InME,,,) = 0, 



(6) 



with the electric conversion K 



Chemical Physiology 

So far in our modeling, the species passing through the 
membrane are H + , HCO^, CI , K + and Na + . Obviously, the 
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Electroneutral 
Cl 



Electrogenic 
2K+* 




Outer medium 



Passive 

Figure 1. Cell model. The cell membrane acts as a capacitance, which is submitted to the membrane potential E,„. Three blocks exchange ions 
between the cytosol and the outer medium: they induce electroneutral transports (with the AE2 and NHE— 1 exchangers), electrogenic currents via 
the Na/K-ATPase, and passive ion channels (for Na + , K + , and Cl~ ). The chemical reactions are assumed to take place within an homogenous 
cytosol. 

doi:1 0.1 371 /journal.pone.0085449.g001 



first two of them are directly involved in the set of protic reactions 
that govern the pH. Since the physiological range of pH lies 
around 7, we must monitor in our analysis the self- ionization of 
water 



Since the dissolved CO2 is very unstable in water, and especially 
in presence of physiological carbonic anhydrase [16], we shall 
directly merge the equilibria (8) and (9) to get rid of C0 2aq and 
obtain an equivalent equilibrium: 



H 2 O^H + + HO - , K n , = [H + ] [HO " 



(7) 



Then, we include the three components of the carbonated 
system. The partial pressure Hco 2 of carbon dioxyde equilibrates 
with aqueous CO2 (according to the Henry law), which undergoes 
two consecutive dissociations. Those reactions and their equilib- 
rium constants are summarized below: 



C0 2gIls ^C0 2aq , 



^=n C o,/[co 2 ] 



(8) 



C0 2glls -H + + HCO3- ,K[ = Ki /K H = [H+ ] [HCO3- ] /Hoo, .(11) 



Finally, we model the protic behavior of all the other species 
within the cytosol by a single equivalent buffer that we name Y, 
acting as 



C0 2aq ^H + +HC0 3 



HCO3 ^H+ +COf' 



^=[H+][HC0 3 -]/[C0 2 ] (9) 



# 2 = [H+][COM/[HC0 3 l (10) 



H H + + Y - , K Y = [H + ] [ Y - ] / [H Y] 



(12) 



with Y 0 = [HY] + [Y]«50-70mM (see [8] for details). Addi- 
tionally, we do not consider the diffusive effects within the cell, by 
assuming homogeneous ionic concentrations. Similarly, the outer 
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media is considered as an infinite bulk with constant physiological 
concentrations of the different entities under consideration. 

Results 

General Theoretical Results 

Self-consistency and modularity. To be self-consistent, our 
model must ensure that each ion and the membrane potential are 
sufficiently maintained by physico-chemical processes (enzymatic 
transport or chemical reaction) in order to avoid some non- 
physical accumulation or discrepancy such as negative concentra- 
tion. Accordingly, we propose a model for a generic cell, restricted 
to the previous components (transporters, ionic permeabilities and 
chemical reactions) which fulfills these criteria. Noticeably, any 
other effector that acts redundantly for pH regulation or is 
expressed in specialized cells can be implemented, provided that 
the above self-consistency is preserved. As an example, we will 
show further how to handle the lactate/H + production and 
transport. The same methodology applies to other mechanisms 
(such as H + -ATPase or Na + -coupled-bicarbonate transporters) 
and any additional chemical reaction. 

Full formal system dynamics. For any reaction indexed by 
k, we will note Xk its associated rate, which is the derivative of its 
chemical extent with respect to time. For each chemical species, 
the concentration temporal derivative is the appropriate summa- 
tion of the chemical molar rates, the exchanger transport rates and 
the passive intake molar rates. The full system is straightforwardly 



{ ce,E m 

d t [K+] 
c?,[Na+] 
0«[CT] 
5,[HO-] 

a«[cof-] 
a,[HY] 

dt[Y~] 
S ( [HC0 3 - 

I e ( [H+] 



•7 r V(/t K + + /. Na + — Acr 
;. K+ +2p NaK 



" Pn&k) 



X 



Na + 



" ^PNaK + PNHE 



Pae 



Xw 
X2 

-Xy 
Xy 

-Pae + x'i-Xi 
-Pnhe+Xw+x'i+Xz+Xy 



(13) 



where /„„ x 2 , and Xy are respectively the molar rate of the 
water ionization (7), the direct formation of HCO-f (11), the 
dissociation of HCO^ (10), and the deprotonation of the 
equivalent buffer (12). The first equation of the above system is 
rewritten from relation (6). 

The main characteristic of protic reactions in water is that they 
have very short relaxation times, from a few microseconds for 
water [1 7] to a few milliseconds for CO2 with carbonic anhydrase 
[16] or without carbonic anhydrase [18]. Since the transmem- 
brane exchanges of protic species (through NHE — 1 or AE2) are 
expected to have characteristic times much larger, we can 
consequendy make the assumption that each protic equilibrium 
is in fact a fast pre-equilibrium. It follows that each involved 
reaction quotient always matches its corresponding thermody- 
namical reaction constant: this can be seen as a set of constraints 
applied to the chemical composition of the aqueous solution. 
Consequendy, if we want to impose a perturbation of this 
composition then those pre-equilibria shall instantaneously pro- 
duce the mandatory chemical extents which ensure that the final 
composition respects the chemical constraints. Accordingly, the 



thermodynamical knowledge of the equilibria (7) to (12) is 
sufficient to solve the kinetic equations (13) in this particular 
biological context. 

In the following, we detail the treatment of the protic reactions 
rates within the equations (13). We show (see Methods in File SI) 
how to derive a set of reduced differential equations for the 5 
dynamic variables E m , [K + ], [Na + ], [CP] and [H + ], within this 
pre-equilibria approximation. 

Steady state characterization. The steady-state values, 
which we note withan asterix, are obtained by setting to zero 
the temporal evolution in the differential system (13) leading to 



r 0 

0 
0 
0 
0 



A 



K + +^Na+ ~^C\- ~PNaK 



^-K+ + 2 ^NaK 
/ ' l Na+ ^^NaK 

PNHE ~~ P*AE ■ 



Pnhe 



(14) 



This system is under-determined since the sum of the first 
equation and of the last two ones minus the second and the third 
one is zero: this is expected since the evolution of E„, is the exact 
conservation of the global electric charge. In our model, the latter 
decomposes into an intrinsic charge Q of all the considered 
components and an excess charge Q xs of all the other "spectator" 
species (proteins, other ions...), leading to CE„, = AQ + AQ XS , 
where A means the difference between the cytosol and the outer 
medium values. It is the integrated form of equation (6). As a 
consequence, the initial condition of the differential system (13) 
gives AQ" 1 . 

In order to determine steady-state values, we first obtain the 
3 



electric equation 0 = 



- A C [- which is here equivalent 



to 



C™ = m 



[cr 



^ 2^K+[ K+ ]+ P Na+[ Na + ]+ P Cl-[Cr 



(15) 



since we only deal with monovalent ions. The relation (15) is the 
Goldman-Hodgkin-Katz potential equation with voltage-depen- 
dent permeabilities and a potential explicitly regulated by Na/K- 
ATPase. 

With I'ae/nhe = ^ae/ ^nhe (see equations (3) and (4)) and 
Lq»1, the last equation of the system (14) can be expressed as a 
function of h* and it reduces to a simple polynomial 



K AE K r fh* 



K[Yico 2 \K r 



( 



\ 



AE/NHE 



1- 



1 



1 ^AE/NHE 



=0(16) 



This analytical relation yields the steady-state pH as a function 
of IIco, an d ^ae/^nhe, since h* is the positive root of equation 
(16). The Figure 2 shows how the exact pH* evolves when those 
parameters are changed and where the acceptable physiological 
limits stand. In particular, for an intracellular pH = 7.2 and 
nco 2 =40mmHg our model predicts Kae/ ^nhe — 0.057: this 
transport ratio matches well the experimental maximal rates of this 
transporters in different systems [19-21]. An interesting feature of 
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our model is the prediction of missing parameters (kinetic and/ or 
thermodynamic) based on the knowledge of steady-state physio- 
logical values (see Results in File SI). For instance, a unique [H + ]* 
is computed from a given Vae, ^nhe an d FIco,- Conversely, the 
ratio Kae/^nhe can be read on Figure 2 from the experimental 
measure of [H + ]* and Ilco,- 

Asymptotic Kinetics Framework 

General philosophy. The set of differential equations (13) 
defines a multiple-scale system (in both time and in concentra- 
tions), since it combines slow chemical rates with fast relaxing 
protic reactions. In either muliple-scale analysis [22] and normal 
forms in central manifolds [23], the slow dynamics are assessed 
around a stationary point. However, in the case used herein, the 
slow dynamics evolve on manifolds generated by the laws of mass 
action (corresponding to each protic pre-equilibrium established at 
its thermodynamical constant) and represent the only valid 
compositions of the system. To the the best of our knowledge, 
this is the first time that a way to compute the constrained 
evolution of the all the involved concentrations has been exposed. 

Chemical system description. If we assume that we have 
N chemical reactions coupling M species X^,., m , then the 
relevant reactions may be written in a generic form, employing the 
braic stoichiometric coefficients V : 



V/el.JV, v,. y X / =0,^, = n /El .. M [X] / v '' 

jeh.M 



(17) 



where K, is the equilibrium constant of the i th reaction. We also 
assume that K, is temporally dependent, so as to reflect the 
possible variations of the external conditions (such as imposed 
changes in partial pressures). 

Fast pre-equilibria consequences. We now assume that 
those N reactions represent fast pre-equilibria. In other words, we 
suppose that the relaxation time of each reaction is infinitely small. 
Accordingly, we define the vector T by its coordinates r, such that 



Vz 



el..N,fi = KiH DC : nj - IT [Xl! iJ = 0. (18) 

V ,-<0 l J V: ,->0 L >J y ' 



Thus, the only permissible evolutions must be satisfied at 
particular time t and for any set of concentrations [X] through the 
relationships: 



f(?,[x]) 
d,f + ^d, [x 



(19) 



where <f> is the Jacobian matrix off with respect to j^xj . 

Response to perturbations. If we perform a small modifi- 
cation <5^xj during St on this system, all the reactions evolve to 

preserve T = 0. The resulting individual chemical extents of 
each reaction form the vector £. We then obtain a modified 
perturbation 



<5[x]=<5[x]+v T ? 
and this in turn must obey 



(20) 



0 = Std t f + <D<5[x ]'. (21) 
The instantaneous chemical extent is readily computed by 



I = - (Ov T ) 1 (dtd t f + <D<5 [x j 



(22) 



We can show that the N x N matrix <t>v T is invertible for any 

admissible set of concentrations. This purely algebraic property 

results from the convexity of the free enthalpies of reactions from 

which the expressions of T and d> are derived, but this purely 

mathematical demonstration is far beyond the scope of this, article. 

Generic asymptotic kinetics. Finally, noting that d, X 

L J slow 

contains global information regarding the "slow"-changing 
variations of all the chemical species, the overall chemical 
evolution of the system is deduced from (22) by setting 



i[x]=<5fa ( [' 



X I . Thus, asymptotic kinetics can be written as 

slow 



<S>8 



L J slo 



8,r) (23) 



thereby illustrating how the fast chemical reactions are damping 
the slow variations. 

Numerical integration. A modular C++ program (available 
upon request) was designed and exactiy encodes the biological 
effectors, the membrane potential and the chemistry into a system 
of algebraically coupled numerical differential equations. The 
integration step was performed by an adaptive Dormand-Prince 
method with a fractional tolerance of 10~ 7 . 

Reduced Model for pH Dynamics 

Evaluation. If it is assumed that all the protic reactions are 
rapid pre-equilibria (see above), then we can derive the proton 
generation rate by using the previous mathematical formalism. 
This formalism provides an algebraic manner to decouple all the 
protic reactions from the catalytic ones explicitly. Implicidy 
however, all protic dynamics may be deduced from the evolution 
of [H + ] which has unfortunately no simple expression. 

Simplified pH dynamics. For a given cell in physiological 
conditions (for which the internal steady pH is around 7.2) we can 
neglect the presence of CO 2- and obtain a slightly simplified rate 
(see Methods in File SI): 



3 ( [H+] «®([H+],Y 0 )( (p AE -Pnhe) + tJJ^H 



0 



[H+] 2 



(24) 



[H+] 2 +ic-„,+^;rico 7 i+/» Y Yo 



with a numerically derived value j8 Y — 6.55 M~ , and for a total 
buffer concentration Yr> Here, p AE and Pnhe depend only on the 
pH and on nco 2 > so that one can simulate the pH with only one 
differential equation. 

Role of the buffer. The steady-state pH* is readily recovered 
from equation (24). The factor 0 emphasizes the preponderant 
role of the chemical couplings pertaining to the evolution of the 
intracellular pH. Indeed, we always have 0< [H + ] 2 /(X[nco 2 ), 
so that for the steady state of a "normal" cell (pH = 7.2, 
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n C02 /40mmHg 




Figure 2. Three dimensional representation of steady-state pH. pH is drawn as a function of the two controlling parameters n C o, and 
Kae/^nhe- The flattest area of this surface stretches over +0.5 units around physiological pH (7.2, in black). This corresponds to values that can be 
reached by cellular systems (red boundaries). 
doi:1 0.1 371 /journal.pone.0085449.g002 



nco 2 =40 mmHg), we estimate ©*<5-10~ 6 . The protons 
dynamics (hence of all the protic species) are sharply damped by 
those chemical couplings. As expected, the buffer causes these 
dynamics to be even further reduced through the term 1 + /SyYo 
in the 0 denominator, up to 30% for a 60 mM buffer 
concentration. 

Natural overshoot. Interestingly, our calculations predict 
that a vanishing physiological protic perturbation will systemati- 
cally produce a pH overshoot around its steady state. Such 
phenomena are well known experimentally [24,25]. The exact 
mathematical demonstration of this phenomenon (see Methods in 
File SI), is valid for acidification or alkalinization, and can be 
applied to model other overshoots observed in different physio- 
logical regulations. 

To explain this in a non-mathematical way, we may perform 
the following thought experiment. Let us assume that a weak 
protonated base HB enters the cell at its steady state. The excess of 
protons produced by the dissociation HB^±H + +B~ is contin- 
uously pumped out of the cell by the regulating enzymes. As a 
consequence, if the cell removes HB from its cytosol, then some 
protons will not neutralize the remaining B . Accordingly, the 
initial pH is reached earlier than expected: the further removal of 
B straightforwardly creates an unexpected depletion of protons 
(basic environment) before returning to the initial situation. This 
describes an overshoot mechanism. 

Steady-state pH: Role of Enzymatic Constants 

We have investigated the changes in the steady-state pH 
resulting from covalent or non-covalent modification of the 
transporters through intracellular signaling cascades, and the 
effects of allosteric activators or inhibitors or mutations that affect 
the transporters parameters. Unless stated, we model these effects 
assuming that a modification is specific and affects only one 



thermodynamical enzymatic constant without changing the others, 
while we keep the kinetic ratio Kae/ ^nhe and rico 2 to their usual 
levels. 

For NHE— 1, it has been shown that within the Monod- 
Wyman-Changeux framework [26], the allosteric constant Lq is 
modified by various stimuli such as growth factor stimulation or 
changes in membrane composition and tension [13,21]. This 
raises the question whether NHE — 1 cooperativity for proton is 
intrinsically important for pH regulation itself and other cellular 
functions. Accordingly, Figure 3 depicts the results of our 
computations of the resulting pH following a modification of Lo- 
Interestingly, it has been shown widely that the activation of 
NHE— 1 by the above-mentioned stimuli [27] decreases Lq by 
one order of magnitude [13], and results in a pH increase of 0.2 to 
0.3 units. Our model yields a pH increase of about 0.3 units for 
Lo, s tim/Lo = 0.\, which is in very good agreement with the actual 
experimental data. We are also able to hypothesize that in some 
cases NHE — 1 might be regulated by altering its microscopic 
affinities for protons. To investigate this, we changed the K r and 
K t constants of NHE— 1 and kept the K,/K t ratio constant, as 
these two affinities correspond to the same site in different 
conformations. Our model predicts that any important change in 
these microscopic affinities would produce a large pH shift, as 
shown on Figure 3. 

For AE2, we investigated either the effect of a variation in the 
affinity K/^ or in the Hill exponent (see Figure 3). We observe a 
less drastic change than those resulting from a NHE — 1 
modification. This makes sense from a physiological point of 
view: due to membrane potential and metabolic activity, cells 
constandy have to compensate for intracellular acidification rather 
than for alkalinization. 

To summarize, we show that pH regulation is robust for two 
main reasons. Firstly, changes in the thermodynamic constants of 
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Figure 3. pH modification by changes in biochemical constants. The reference values are L 0 = 10 3 , ST r = 1.810 -8 M, AT, = 3.6 10 6 M, 

Kxe = 10 mM and y= 1. For K, and K, modification, the ratio K r /K, is held constant (see text for the explanation). 
doi:1 0.1 371 /journal.pone.0085449.g003 



pH regulation systems, that in vivo could arise from mutations or 
from interspecific variations, induce minute modifications of the 
steady-state pH. Secondly, pH can relax back to its physiological 
value, because changes in constants are very easily overcome by 
slight modifications of the maximal rates of the transporters, i.e. the 
amount of transporters expressed at the plasma membrane. This a 
posteriori validate the hypotheses and choices described in the 
introduction. 

Characteristic Time Scales in pH Regulation 

In order to investigate the dynamics, we integrated the 
differential system (13) as described previously. We choose to 
approximate the average fibroblaste shape using a prolate 
spheroid model of length 25^m and a diameter of 10 /mi, leading 
to a surface 5 = 653 /an 2 and a volume V= 1309 /an 3 . 

We found an anionic charge excess of AQ XS ~ — 30 mM for this 
configuration, which mainly corresponds (i) to the excess of 
negative charges found on the surface of intracellular proteins and 
(ii) to the bulk of negative charges provided by the first dissociation 
of phosphate groups [8] . 

Relaxation times around the steady-state values. We 
performed the linear stability analysis of the differential system (13) 
which also provides the relaxation constants of the independent 
variables. We deduced the raw and typical relaxation time 
constants for our cell model by setting the equivalent buffer 
concentration to zero. Firstly, we obtain a 3 ms characteristic time 
that predominantly corresponds to the relaxation of E m : obviously 
the membrane potential adjusts itself very quickly to a change in 
the ionic composition, but since it does not produce chemical 
species per se, it does not influence the chemical rates. 

Secondly, we have two similar time constants representing the 
relaxation of a perturbation of all the ions within 8 and 15 
minutes. 



Finally, for a fixed II^, the perturbed concentration <5[H + ] 
dynamics obeys to 



with 1 jxh - 
©*«510- 



0 



d,<5[H+] = 
<5(Pnhe"Pae) 



-5[H+]/% h 



(25) 



3[H + ] 

the proton relaxation time T/, 



Since the relation (24) provides 



reaches about 5 



minutes, which is consistent with the experimental observations 
from a plethora or reports [24,28,29] . 

Illustration on Pathophysiological Situations 

Forced acidosis: IT.CO2 and intracellular buffer. We 

simulated an artificial increase of Tlco 2 of 20% during 1 minute 
(from 40 mmHg to 48 mmHg) followed by a return to the normal 
within 5 minutes. The resulting pH dynamics with and without the 
equivalent buffer HY/Y~ are presented in Figure 4A. As shown 
in Figure 4B, the HCO^ excess increases the CI intake, and the 
acidification increases the Na + intake, while the K + level is 
remarkly stable as expected. The net ionic currents produce a 
concurrent tiny 2 fiV depolarization. We note that, via the 
chemical couplings, the aprotic species dynamics are also 
dampened by the presence of the buffer. As expected, the different 
timescales are also respected once the perturbation is over and all 
the concentrations are relaxing towards their steady-state value: 
the pH needs only a few minutes to recover its physiological level, 
while the other ions rather require tens of minutes to reach their 
final balance. The numerical simulation also points out the 
predicted overshoot with and without buffering. 

If IIco2 is held to at its maximum increase, then the pH 
converges to a new value that can be deduced from Figure 2. In 
such a case, the pH curve is similar to Figure 4A except that it 
converges monotonously towards almost 7.18 after the initial 
decrease, and no overshoot occurs. At the same time, the other 
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Figure 4. Forced acidosis by a simulated IIco 2 spike. The rise and the decrease of n C o, are highlighted by blue areas. (A) The expected pH 
overshoot takes place with or without the presence of a 60 mM physiological buffer (dashed lines). (B) The ionic ratios relative to the initial values are 
reported as well. 

doi:10.1371/journal.pone.0085449.g004 



ions find a different balance. The results of this specific fIco 2 j um p 
are shown on Figure S2A&B in File SI. 

High-flow lactic ischemia. Here we show how to expand 
the model in order to probe the consequences of a slight hypoxia 
without ATP depletion. We model it with a lactic acid production, 
while we hold the enzymatic constants and flccb to their normal 
values. Since our model is modular, we first consider the 
dissociation of the lactic acid LaH, namely LaH?±La~ +H + 
with K a = 10-™ 6 . The cell removes lactates and their accompa- 
nying protons (1:1) through monocarboxylate transporters (MCT) 
[30,31]. The latters follow a Michaelian law defined by 
K m ~ 30 mM and by an observed maximum rate 
a_ ~ 1 mM/min [32]. The global lactic acid production is around 
a + ~ 1 mM/min in an hypoxic skeletal muscle [33]. Consequent- 
ly we simply have (i) to append 



3, [LaH] 
5 t [La-] 



-Xa + Pa 
+ Xa 



(26) 



where % a * s me molar rate of the lactic acid dissociation and (ii) to 
include y ul to 3j[H + ] within the differential system (13). Here we 
impose 



Pu- 



if 0</<r fl 



-<x_[La-]/(£» + [La-]) iff>J fl 



(27) 



where T a is the ischemia duration. The resulting pH is shown in 
Figure 5A. 

Consequently, an extra term appears in the overall protonic 
rate: 
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Figure 5. Simulated ischemia. Lactate production occurs with a constant Ilco,/ a 60 mM buffer, and with transporters working at their nominal 
level. The lactate accumulation is highlighted by a green area. (A) The predicted pH overshoot takes place during the lactate removal by the 
monocarboxylate transporters. (B) The corresponding ionic ratios relative to the initial values. 
doi:1 0.1 371 /journal.pone.0085449.g005 



?,[H+]cr(^, + [H+])(p NHE -p AE )+^,p u . (28) 

The acidification of the cytosol produces a massive Na + 
overload (see Figure 5B), which is experimentally observed [28] 
and corresponds to a stimulated action of NHE — 1 . For a fixed 
nco 2 , die fall of [HCO^] during the lactate production decreases 
the chloride intake, as shown in Figure 5B. The net ionic currents 
induce a small 200 /iV hyperpolarization over the simulation. As 
demonstrated in the Methods of File S 1 , the structure of equation 
(28) leads to an expected pH overshoot, which occurs after 
t> T a = 15 minutes. The fast regulating couple formed by 
NHE — 1 and AE2 allows the pH to follow the rate limiting 
lactate expulsion. 



Discussion 

Adequacy with Experimental Data 

As previously stated, the main purpose of this study was to build 
a mathematical depiction of intracellular pH regulation and 
investigate whether it had analytical solutions and produced 
biologically relevant simulations. This last section intends to 
further challenge our study by confronting real experimental data. 
To avoid potential biases, we decided against the use of our own 
data and instead to choose one of the pioneer experiments within 
the large body of published intracellular pH measurements 
generated by independent groups in the last four decades. 
Namely, we use here experimental recordings performed in one 
of the chief studies on intracellular pH regulation published by 
Roos and Boron in 1981 [24]. In this study (figure 5 A of the 
original article), a Helix neuron was submitted to a 10 minutes 5% 
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CO2 pulse. Its pH dropped from 7.35 to 6.85 and returned to 
normal after the pulse, with a noticeable overshoot. Details of the 
calculation and graphical results produced by our simulation are 
given in Materials and Figure S3A in File SI. Taken together, they 
show that only very minimal modifications of the constants of the 
system, well within differences found between different cell lines 
such as fibroblasts and neurons, have to be applied to converge to 
the resting intracellular pH measured in experimental conditions 
and that very satisfactory matches are obtained between calculated 
and experimental values for nco 2 an d intracellular pH. 

Main Outcome 

This study describes the first fully coupled and self-consistent 
mathematical system for intracellular pH regulation. For this, we 
constructed a minimal system that is uniquely based on the kinetic, 
electric and chemical equations describing the molecular processes 
pertaining to intracellular pH. This strategy is very different from 
classical heuristic methods used to model biological processes, that 
are mostly built on phenomenological equations deduced from 
fitted curves. It also avoids the introduction of ad hoc fluxes and/ or 
constants to ensure the convergence of the numerical simulations 
with experimental data. Importandy, the present approach allows 
analytical processing. It shows, for the first time, that the dynamics 
of pH can be described by a reduced differential equation, and 
that steady-state intracellular pH values are in fact analytical 
solutions. Besides, despite the formal complexity provided by the 
large body of equations used here, the calculated numerical values 
of pH, ionic concentrations and membrane potential converge 
towards physiological values, with time evolutions that are very 
reminiscent of experimental behaviors. The last remarkable 
finding is the demonstration that any additional phenomenon 
that directly or indirectly impacts pH can be mathematically 
included without violating our model, provided that its equation is 
not ill defined. At this step, it is important to notice that here, we 
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